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Abstract 

We discuss kinematic methods for determining the masses of the particles in events at a hadron 
cohider in which a pair of identical particles is produced with each decaying via a series of on- 
shell intermediate beyond-the-SM (BSM) particles to visible SM particles and an invisible particle 
(schematically, pp — )• ZZ + jets with Z — > Aa — t- Bba — )• Ccba —)•...—)• cba . . . + N where 
a,b,c, . . . are visible SM particles or groups of SM particles. A, B,C, . . . are on-shell BSM particles 
and is invisible). This topology arises in many models including SUSY processes such as squark 
and gluino pair production and decay. We present the detailed procedure for the case of Z — > 
3 visible particles + N and demonstrate that the masses obtained from the kinematic procedure 
are independent of the model by comparing SUSY to UED. 
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I. INTRODUCTION 



Many solutions to the hierarchy problem require new particles whose loop corrections 
to the Higgs mass-squared cancel the quadratically divergent Standard Model (SM) loop 
corrections. The masses of the new particles (especially the particle or particles which 
cancel the SM top loop) should be sufficiently small that the Higgs mass can be naturally 
below the TeV scale. On the other hand, if the new particles have masses below 0{TeV) 
then LEP observables will be strongly affected if they are exchanged at tree level or can 
be singly produced. Both these latter possibilities are automatically removed if there is 
a symmetry under which the new particles are odd and the SM particles are even. In 
particular, the new particles can then only contribute to the electroweak observables at the 
loop level, and new particles with masses of order a few hundreds of GeV can be compatible 
with the data. In such scenarios the lightest of the new particles is automatically stable and 
it should be neutral for consistency with bounds on new charged stable matter. Typically, 
it is also weakly interacting. Such a weakly interaction massive stable particle (WIMP) will 
be invisible, leading to "missing energy" in particle detectors. This scenario is also highly 
desirable since such a WIMP can readily provide the dark matter known to be present in 
the universe. 

Almost all the models with dark matter candidates also contain additional particles that 
are charged not only under the new symmetry but also carry SM "charges," most often 
including color. At a collider, these new particles must (and will) be pair-produced, and 
since they are heavier than the dark matter particle, they will cascade decay down to it. In 
many cases, this cascade radiates SM particles in a series of A — > Be, 1 — body — > 2 — body 
decays, in which A and B are new physics particles while c is a SM particle. (In some 
cases, phase space restrictions force one of the new particles off-shell and A ^ B*c ^ Cdc, 
1 — body — )■ 3 — body decays are relevant.) Since the final step in the chain will yield a dark 
matter particle, the typical collider signals for such a scenario will be jets and/or leptons 
plus missing energy. 

Supersymmetry (SUSY) is the most popular model of this type. In SUSY, the new 
symmetry is termed matter Parity (sometimes called i?-parity). Its conservation implies that 

the Lightest Supersymmetric Particle (LSP) is stable. In most supersymmctric models the 
LSP is the lightest neutralino, which is a good dark matter candidate. It appears at the end 
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of every supersymmetric particle decay chain and escapes the detector. All supersymmetric 
particles are produced in pairs, resulting in at least two missing particles in each event. 

Other theories of TeV scale physics with dark matter candidates have been recently 
proposed. They have experimental signatures very similar to SUSY: i.e. multi leptons 
and/or jets plus missing energy. For instance, Universal Extra Dimensions (UEDs) 

little Higgs theories with T-parity (LHT) [sj], and warped extra dimensions with a 
parity j4| belong to this category of models. 

Clearly, being able to reconstruct events with missing energy is an important first step 
to distinguish various scenarios and establish the underlying theory. In addition, studies |5| 
suggest that the mass of the dark matter particle, and the masses of any other particles 
with which it can coannihilate, need to be determined to within a few GeV in order to be 
able to compute the dark matter density in the context of a given model. A very important 
question is then whether or not the LHC can achieve such accuracy or will it be necessary 
to wait for threshold scan data from the ILC. The goal of this paper will be to provide 
details regarding the kinematic techniques developed in Refs. js, 7| that provide the needed 
accuracy using just LHC data. For the case of 3 visible particles per decay chain, the focus 
of this paper, we also show that the kinematic technique gives masses for the BSM particles 
that are completely insensitive to the particular model by comparing a SUSY case to a UED 
case where the decaying BSM and final invisible BSM particles in the two cases have the 
same masses. This implies that it is unnecessary to determine the overall mass scale of the 
BSM particles using model-dependent information, such as total cross sections. Indeed, to 
fully test a potential model, it is necessary to first determine the masses of the produced 
particles just based on kinematic information. Once the masses are known, there are many 
chain decay configurations for which it will be possible to determine the four-momenta of 
all the particles on an event-by-event basis. The four-momenta can then be employed in 
computing the matrix element squared for different possible spin assignments. In this way, 
a spin determination may be possible which, in combination with cross section information, 
can be used to distinguish different models. 

In recent years there have been numerous studies in the context of SUSY-like theories 
of how to measure the super-partner masses just based on kinematic information 
In some cases the procedures employ a single long decay chain of super-particles, usually 
requiring 3 or more visible particles in the decay chain in order to have enough invariant 
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FIG. 1: A decay chain in SUSY. 

mass combinations of the visible particles to achieve sensitivity to the absolute mass scale, 
as opposed to simply mass differences. Squark decay, Fig. [H is an example of one such chain. 

Our approach has been to pursue alternative procedures that employ events in which a 
pair of identical (particle and antiparticle, e.g. squark plus anti-squark) BSM particles is 
produced and both decay in the same manner. In such an event, information from both 



decay chains in the event can be included at once. In our first paper jo], we tackled the 
difficult case where we assumed that only two particles appeared in each chain decay, e.g. 
making use only of the leptons appearing in Fig. [H In this case, kinematic constraints alone 
can not give a discrete solution for the unknown masses. Nonetheless, the space of the 
allowed solutions does contain enough information about the new particle masses and they 



can be extracted using a statistical procedure iGl. The mass determination can be further 

nn 

improved by combining with other kinematic variables such as Mt2 |9|, |23[ . However, these 
kind of analyzes usually require large statistics in order to achieve a reasonable precision. In 
this paper, we provide details on the case where 3 visible particles are present in each chain 
decay. In this case, a single pair of events provides enough information to yield a discrete 
set of possible masses. Therefore, very few events are needed for a rough determination of 
the masses, in contrast to the statistical methods which rely on the availability of a large 
number of events. 

The general topology on which we focus is then that of Fig. [2l After including com- 
binatorics and resolution, we will achieve root-mean-square (rms) accuracies on the three 
underlying masses in the decay chain of order a few GeV (depending upon the number of 
available events) with a systematic shift that can be easily corrected for. This result is fairly 
stable when backgrounds are included so long as S/B > 2. 

The organization of the paper is as follows. In Sec. HIl we give the general counting of 
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FIG. 2: The event topology we consider. 

constraints and unknowns for single chain and multiple chain events. In Sec. IIIIl we give 
a more detailed exposition regarding solving the topology of Fig. |2l In Sec. IIVI we first 
demonstrate how the masses of the Z, Y, X and particles in Fig. |2]can be very precisely 
determined using just a few events if there are no effects associated with combinatorics, 
particle momentum measurement resolutions or backgrounds. We then develop the very 
crucial strategies for dealing with the realistic situation where combinatorics, resolution 
effects and backgrounds are present. We still find good accuracies for all the masses using 
only the kinematic information contained in the available events. We study the accuracy of 
the mass determinations as a function of the available number of events and as a function of 
the signal to background ratio. In sec. |Vl we compare results for the SUSY and UED cases 
and show that the masses determined are independent (to within one to two GeV) of which 
model is employed. We summarize and present additional discussion in Sec. IVII Some of 
the material in sec. [inland sec. IIV]has appeared in Ref. |7l], but is included in the present 
article for completeness and to simplify some of the discussions. 

II. CONSTRAINTS COUNTING 

To begin, it is useful to perform a general counting of observables and constraints for 
various different configurations. We consider first the counting when only one decay chain 
in the event is considered at a time. We then show the increase in constraints possible if 
both decay chains in each event are considered at once. 
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A. Single decay chain case 



We begin with the chain decay X — aY — )■ abN, where the 4-momenta of the SM- 
particles a and b are directly measured. For each event, there are four unknowns due to 
the unobserved 4-momentum of the A^. In addition, we have the three unknown masses 
(the same for every event) mx, ttly, and rriN. These are subject to 3n constraints coming 
from requiring that the X, Y, and be on their mass-shell. Thus, after n events we have 
3 + 4n — 3n = 3 + n free parameters. No matter how many events we examine, we will not 
be able to obtain a discrete solution (or set of solutions) for nix, ^y, and rriN. 

Next, consider a decay chain with three observable SM particles: Z — t- aY abX 
abcN. In this case, the number of unknowns after n events is 4 + 4n and the number of 
constraints is 4n (the Z, Y, X and masses, which are the same for every event). After n 
events we then have 4 + 4n — 4n = 4 free parameters, which basically correspond to the 4 
unknown masses of the decaying particles. In this case, each event will determine a region in 
the {mz,mY,mx,mN) mass space and, as more events are accumulated, in an ideal world 
the region of mass space consistent with all events would become more and more restricted, 
but it will never reach a discrete point (or a set of discrete points) with any number of events. 
To pin down the actual mass point, one needs to use additional information by examining 
the end points of certain kinematic distributions. 

If one considers (as in [3]) a chain with four fully measured SM-particles, A — )■ aZ — )■ 
abY — )■ abcX — )■ abcdN, then we have 5 on-shell particles, and after n events, we end up 
with 5 + 4n unknowns and 5n constraints. The number of free parameters is then 5 — n, im- 
plying that (up to discrete ambiguities associated with a high order polynomial and ignoring 
combinatorics and resolution) n = 5 events would be sufficient to solve for the resonance 
masses. However, combinatorics and resolution effects will considerably complicate the sit- 



uation, as already apparent from the study of [13|, where they assume that my, mx and 
rriN are known, leaving, in principle, 2 + 4n unknowns and 5n constraints after n events, 
implying that only n = 2 events would be needed to solve for the remaining 2 masses, 
and mz- After including combinatorics and resolution, Ref. 13(| needed many more than 2 
events in order to get a mass determination. 

The general counting procedure is apparent. If the decay chain has A^a — 1 on-shell 
decaying particles and a final invisible particle, then after n events there will be 4n unknowns 
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associated with the 4-momenta of the invisible particle. There will also be Na unknowns 
corresponding to the unknown masses of all the BSM particles. These unknowns will be 
subject to Nj^n constraints from the requirement that all the BSM particles have the same 
on-shell masses in each event. The number of unknowns after n signal events will then be 



For A^^ < 4, a discrete solution or set of solutions is not possible regardless of how many 
events are available. The actual masses may still be obtained by combining additional 
information from the kinematic distributions for A^^ = 4 {e.g. the squark decay case of 
Fig. [1]). For Na = 5 {e.g. the SUSY g ^ qq ^ 1(1X2 ~^ 11^^ ~^ QQ^^Xi decay chain) n = 5 
events will give a discrete set of solutions for the masses. Still longer decay chains would 
require fewer events to obtain a set of discrete possibilities for the BSM particle masses. An 
example of a longer decay chain with A^^ = 6 (implying that 3 events would give a discrete 
set of solutions for the masses) would be 



where the Z would be observed in one of its visible decay modes. By considering the full 
event at once through inclusion of information from both decay chains, one need not resort 
to such long decay chains in order to get to the point of having a discrete set of solutions 
for the masses. 

B. Using the whole event, i.e. both decay chains 

When considering the whole event at once, the constraint counting proceeds differently. 
Assuming that there are two invisible particles present in the final state, the number of 
unknowns associated with their 4-momenta in all n events is 8n. Requiring that the sum 
of the invisible transverse momenta equal minus the sum of the visible transverse momenta 
imposes 2n constraints on these unknowns. In addition, let us suppose that the topology 
is such that Nb masses are unknown {Nb includes the unknown masses of the invisible 
particles, which we do not require to be the same at this point in our counting). Each event 
will also be subject to a number A^^ of on-shell mass constraints, including the requirement 
that the two invisible particles have masses equal to their on-shell (unknown) values. Then, 



Nu = NA + ^n- Nati = A^^ + (4 - A^^)n . 



(1) 



g ^ qq ^ qqXs ^ qqZx2 ^ QQ^li qqZUxi 



(2) 
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after n events there are Nau constraints. Thus, the number of unknowns before imposing 
constraints is Nb + Qn and the number of constraints is Nau, leaving 



Nu = Nb + {6- NA)n (3) 

unknowns. If we consider an event with A^^ — 2 > 4 on-shell decays and require that there 
be no unknowns, i.e. Nu < 0, after ns events, we find 

ns>^^. (4) 

Of course, in general Nb < Na. For symmetric chains, Nb = N aI'^- For chains in which 
only the final missing particles are assumed to have the same mass, Nb = Na — 1. The 
particular case we focus on in this paper. Fig. O corresponds Na = 8 and Nb = 4. In 
this case, ns = 2 events will lead to Nu = 0, implying a discrete set of solutions for the 
unknown masses of Z, Y, X, N. Were we to consider the case of two identical chains with 
only 2 visible particles in each and unknown masses for Y, X, N, one would have A^^ = 6 
and Nb = 3, leading to Njj = 3 (corresponding to the unknown Y,X,N masses). This 
is the case considered in [6|; a discrete set of solutions will never emerge with any finite 
number of signal events. Additional information from kinematic distributions is needed to 
pin down the masses. Conversely, if one goes to symmetric longer decay chains (such as 
that of Eq. ([2])) with A^^ = 12 and A''^ = 6, then just one event will be sufficient to give a 
discrete set of mass solutions. The general problem with longer decay chains is that they 
are less likely to occur and harder to identify even if they occur, as the visible particles may 
be lost or hard to isolate. The shortest decay chains that can give rise to discrete solutions, 
i.e., the topology of Fig. [21 with A^^ = 8 and A''^ = 4 yielding ns = 2, is likely to be the 
optimal configuration for the mass determination. 

Extension to more missing particles, or quadratic constraints is straightforward, but 
increases the order of the equations to be solved, requiring more advanced polynomial solvers. 
It also requires that more masses be specified. For example, let us assume some definite 
topology for which each event contains 3 'missing' particles. If we do not presume any 
mass equalities among the total number, Na, of decaying resonances and missing particles, 
then the number of unknowns after n events is Na + (4x3 — 2)n, where the —2 is the 
transverse momentum constraint setting the visible transverse momentum equal to minus 
the invisible transverse momentum in each event and the 4 x 3 = 12 just corresponds to 
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the 4 unknown components of each of the 3 invisible particles' 4-momenta. The number of 
mass-shell constraints is (by definition) Nau. If the topology is such that Na > 10 on-shell 
masses can be reconstructed from the visible momenta and the invisible momenta and if 
only Nb of the A^^ masses are independent, then the number of unknowns after n events is 
Nb + (4x3 — 2)n and the number of constraints is Nau. Again neglecting possible relations 
among these masses and requiring the final number of unknowns, Nu = Nb + lOn — Nau, 
to be or negative in order to have fewer unknowns than constraints after n = ns events, 
one finds that 

events would lead to a certain set of discrete mass solutions. If all the invisible particles are 
the same then Nb < Na — 2. 

We have classified many possible decay chains which fall into a category such that a 
small handful of events could give an essentially unique mass spectrum in the absence of 
combinatorial and experimental resolution effects. However, generically speaking one wishes 
to keep the chains as short as possible while consistent with a small number of events being 
sufficient to yield a discrete spectrum of mass solutions. This is because (a) shorter chains 
are easier to isolate on an event-by-event basis and (b) combinatorial and resolution smearing 
of the solutions may be lessened. 



III. BASIC EQUATIONS FOR THE TOPOLOGY OF FIG. 2 

The topology on which we focus in this paper is that given in Fig. |2l As sketched in 
the previous section, this is an ideal topology for precise mass reconstruction. Assuming 
TTiN = Tn^'i mx = Tnx'i rny = my, mz = mz' and denoting the 4-momenta for particles 
i{i = 1 . . . 8) with Pi, we have 

pl=pl(=m%), (6) 
iPl+P3f = iP2+P4Yi=mj,), (7) 
(Pi +P3+ Pbf = iP2 +P4+ P&f{.= m^), (8) 

{Pl + P3 + P5 + Plf = iP2 +P4+P6+ P^f{,= ml). (9) 
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We assume further that the only invisible particles are particles 1 and 2, and thus have two 
more constraints, 

P1+P2 = Pmiss, P\+PI= Plaiss- (10) 

There are 8 unknowns in Eqs. ([6]) through (fl0|) . namely, the 4-momenta pi and p2 of the 
missing particles. Therefore the system is underconstrained and we cannot solve the equa- 
tions. This situation changes if we add a second event with the same decay chains. Denoting 
the 4-momenta in the second events as = 1 ... 8), we have 8 more unknowns, qi and ^2, 



but 10 more equations, 

ql = ql=pl (11) 

(gi + qz? = (92 + q^y = iP2 + Pa)\ (12) 

{qi + 53 + q^f = iq2 + q4 + q&f = iP2 +P4 + p&f, (13) 

(gi + gs + + qif = {q2 + g4 + ge + qsf = (p2 +P4 + P6 + psY, (14) 

gi + ^2 = qmissj qf + ql = qmiss- (i5) 



Altogether, we have 16 unknowns and 16 equations. The system can be solved numerically 
and we obtain discrete solutions for pi, p2, gi, and q2 and thus the masses m^, mx, my, 
niz- Note that the equations always have 8 complex solutions, but we will keep only the real 
and positive-energy ones which we simply call "solutions" in the rest of the paper. Thus, 
up to a certain number of discrete ambiguities we can determine the Z, F, X, N masses by 
pairing any two signal events. Even a few pairs of events are typically sufficient to eliminate 
the discrete ambiguities due to higher order equations. However, effects such as wrong 
combinations and solutions, initial and final state radiation, experimental resolutions, and 
background events will add complications, which we address in Sec. IIVI 

The equations ([6]) through f|T5|) can be easily reduced to 3 quadratic equations plus 13 
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linear equations, 

2 2 2 2 /ia\ 
Pl=P2 = Ql=Q2^ (16) 

2pi -^3+^3 = 2P2 ■P4+pI = 2gi ■ g3 + ^3 = 2^2 ■ 94 + ql, (17) 
2(pi + Ps) ■P5+pI = 2(p2 + P4) ■P6+pI = 

= 2(gi + 93) ■ gs + gs = 2(92 + 94) ■ ge + gl (18) 

2(pi + P3 + Pb) ■P7+P7 = 2(p2 + P4 + Pe) ■P8+pI = 

= 2(gi + g3 + gs) ■ g? + g? = 2(g2 + g4 + ge) ■ gs + ql, (19) 

P?+P2=Pmiss, P\+PI=PLss, (20) 

gr + g2=gmi.., gi+g2=gL.., (21) 

where all but the first line are linear equations because ^3,4,5,6, 7,8 and g3,4,5,6,7,8 are all visible 
measured momenta. In general, the above equation system has 8 complex solutions, each of 



which could be real. This can be shown by calculating the Grobner basis [28|, in which the 
system is transformed to an 8th order univariate equation plus 15 linear equations. Since 
the other 15 equations are linear, it is straightforward to solve for the other 15 variables 
once the 8th order equation is solved. Commercial software such as Mathematica uses this 
method. However, it consumes an intolerably long time for a single or small number of 
PCs. We take a simpler and faster approach which is described in detail in the Appendices. 
In our method, instead of ending up with an 8th order equation, we obtain a 9th order 
univariate polynomial equation and therefore introduce a fake solution in addition to the 
true solutions. The 9th order univariate polynomial equation is numerically solved using 



algorithm TOMS/493 311] • The fake solution can be easily eliminated by substituting back 



all solutions in the original equations. 



IV. APPLICATIONS 



A. SUSY point SPSla 



For illustration and easy comparison to the literature, we apply our method for the SUSY 



point, SPSla 32], although many of the discussions below apply for generic cases. For 
SPSla, the particles corresponding to N,X,Y,Z are x?, ini^i = e//x), Xg, gL(g = d,u,s,c) 
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FIG. 3: We plot the number of mass solutions (in 1 GeV bins — the same binning is used for the 
other plots) vs. mass in the ideal case. All possible pairs for 100 events are included. Signal events 
only. 

respectively. The masses are 



with the final two numbers corresponding to up/down type squarks respectively. Since 
rrir 'T^p./Ti the i = t case is an important background. We generate events with PYTHIA 



We first consider the ideal case: no background events, all visible momenta measured ex- 
actly, all intermediate particles on-shell and each visible particle associated with the correct 
decay chain and position in the decay chain. We also restrict the squarks to be up-type 
only. In this case, we can solve for the masses exactly by pairing any two events. The only 
complication comes from there being 8 complex solutions for the system of equations, of 
which more than one can be real and positive. Of course, the wrong solutions are different 
from pair to pair, but the correct solution is common. The mass distributions for the ideal 
case with 100 events (no kinematic cuts applied) are shown in Fig. [31 Note the logarith- 
mic scale. As expected, we observe 5-function-like mass peaks on top of small backgrounds 
coming from wrong solutions. On average, there are about 2 solutions per pair of events. 

The (5-functions in the mass distributions arise only when exactly correct momenta are in- 
put into the equations we solve. To be experimentally realistic, we now include the following. 



TUN = 97 A GeV, mx = 142.5 GeV, 

my = 180.3 GeV, mz = 564.8/570.8 GeV , 



(22) 



6.4 |35|. 
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mass (GeV) 

FIG. 4: Number of mass solutions versus mass after including all combination pairings for 100 
events. Signal events only, with only combinatoric ambiguities included. 

1. Wrong combinations. For a given event a "combination" is a particular assignment 
of the jets and leptons to the external legs of Fig. O For each event, there is only one correct 
combination (excluding 1357 ^ 2468 symmetry). Assuming that we can identify the two 
jets that correspond to the two quarks, we have 8 (16) possible combinations for the 2//2e 
(4yU or 4e) channel. The total number of combinations for a pair of events is the product 
of the two, i.e. 64, 128 or 256. Adding the wrong combination pairings for the ideal case 
yields the mass distributions of Fig. HJ Compared to Fig. El there are 16 times more (wrong) 
solutions, but the 5-function-like mass peaks remain evident. 

2. Finite widths. For SPSla, the widths of the intermediate particles are roughly 
5 GeV, 20 MeV and 200 MeV for q^, X2 ^^"^ ^R- Thus, the widths are quite small in 
comparison to the corresponding masses. 

3. Mass splitting between flavors. The masses for up and down type squarks have a 
small difference of 6 GeV. Since it is impossible to determine flavors for the light jets, the 
mass determined should be viewed as the average value of the two squarks (weighted by the 
parton distribution functions). 

4. Initial/final state radiation. These two types of radiation not only smear the 
visible particles' momenta, but also provide a source for extra jets in the events. We will 
apply a pr cut to get rid of soft jets. 

5. Extra hard particles in the signal events. In SPSla, many of the squarks come 
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from gluino decay (g — )■ qqi), which yields another hard q in the event. Fortunately, for 
SPSla rrig — m-qj^ = 40 GeV is much smaller than 7715-^ — m^o = 380 GeV. Therefore, the q 
from squark decay is usually much more energetic than the q from g decay. We select the 
two jets with highest px in each event after cuts. Experimentally one would want to justify 
this choice by examining the jet multiplicity to ensure that this analysis is dominated by 
2-jet events, and not 3 or 4 jet events. 

6. Background events. The SM backgrounds are negligible for this signal in SPSla. 
There are a few significant backgrounds from other SUSY processes: 

(a) (Jl — )■ qX2 ~^ 1'^^ ~^ qT'TXi ^oi one or both decay chains, with all r's decaying 
leptonically. Indeed, X2 ~^ '^^ largest partial width, being 14 times that of xH ~^ N^- 
However, to be included in our selection the two r's in one decay chain must both decay to 
leptons with the same flavor, which reduces the ratio. A cut on lepton px also helps to reduce 
this background, since leptons from r decays are softer. Experimentally one should perform 
a separate search for hadronically decaying tau's or non-identical-fiavor lepton decay chains 
to explicitly measure this background. 

(b) Processes containing a pair of sbottoms, which have different masses from the first 
two generations. Since h jets are distinguishable, a separate analysis should be performed 
to determine the h squark masses. However, this presents a background to the light squark 
search since 6-tagging efficiency is only about 50% at high pt- 

(c) Processes that contain a pair of x^'s, not both coming from squark decays. For these 
events to fake signal events, extra jets need to come from initial and/or final state radiation 
or other particle decays. For example, direct x^ pair production or X2 + 5' production. These 
are electroweak processes, but, since xH has a much smaller mass than squarks, the cross- 
section is not negligible. In our SPSla analysis, the large jet pt cut reduces this kind of 
background due to the small rrig — rriq^ . 

7. Experimental resolutions. In order to estimate this experimental effect at the LHC, 
events in both signal and the aforementioned SUSY backgrounds are further processed with 



PGS 1381]. Note that in [7], we used ATLFAST for the detector simulation. Compared with 
ATLFAST, PGS has more stringent lepton isolation cuts, therefore we obtain fewer events. 
Nevertheless, as shown below, the results turn out to be similar. All objects including jets, 
isolated leptons and missing px are taken directly from PGS. 
The cuts used to isolate the signal are: 




n 
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FIG. 5: Mass solutions with all effects 1-7 included and after cuts I - 111 for the SPSla SUSY 
model and L = 300 fb"^ Ah effects incorporated, including backgrounds. 

I) 4 isolated leptons with pt > 10 GeV, |?7| < 2.5 and matching flavors and charges 
consistent with our assumed X2 ~^ ^ ^ Xi decay; 

II) No 6-jets and > 2 jets with px > 100 GeV, \ri\ < 2.5. The 2 highest-pT jets are taken 
to be particles 7 and 8; 

III) Missing pr > 50 GeV. 

For a data sample with 300 fb~^ integrated luminosity, there are about 620 events left after 
the above cuts, out of which about 420 are signal events. After taking all possible pairs for 
all possible combinations and solving for the masses, we obtain the mass distributions in 



From Fig. [5], we see that the mass peaks are smeared but still present around the input 
masses. The analytical formula for the distributions are unknown, so we estimate the masses 
by reading the peak positions. To minimize the effect from statistical fluctuations, we fit 
each distribution using a sum of a Gaussian plus a (single) quadratic polynomial and taking 
the maximum positions of the fitted peaks as the estimated masses. We will use this function 
as the "standard fit" throughout this article. The fitted range is restricted to be above the 
half height. The fitted curves are superimposed on the mass distributions in Fig. [5l which 
yields {78.4, 134.2, 181.5, 553.9} GeV for the masses. Averaging over 20 different data 
samples, we find 



Fig. El 



= 76.7 ±2.0 GeV, = 134.6 ± 2.2 GeV, 

my = 178.9 ± 3.8 GeV, mz = 561.6 ± 5.4 GeV. 



(23) 
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FIG. 6: For each event, and each combination, c, associated with that event, we count the 
number, Npair{c,i), of events that can pair with it and give at least one solution. The plot shows 
the frequency of occurrence of different values of Npairic,i). All effects are incorporated, including 
backgrounds. The plot is for the SPSla case (for which the total number of signal+background 
events is 620 for L = 300 fb"^). In the bias reduction procedure, any choices of c, i yielding 
Npair{c,i) to the left of the red line (corresponding to 75% of the total number of events) are 
discarded. 



The statistical uncertainties are very small, but there exist biases, especially for the two 
light masses. In practice, we can always correct the biases by comparing real data with 
Monte Carlo. Nevertheless, we would Uke to reduce the biases as much as possible using 
data only. In some cases, the biases can be very large and it is essential to reduce them 
before comparing with Monte Carlo-we will see an example later. 

The combinatorial background is an especially important source of bias since it yields 
peaked mass distributions that are not symmetrically distributed around the true masses, 
as can be seen from Fig. |H This will introduce biases that survive even after smearing. 
Therefore, we concentrate on reducing wrong solutions. 

First, we reduce the number of wrong combinations by the following procedure. For each 
combination choice, c, for a given event, i {i = l,Nevt), we count the number, Npair{c,i), 
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of events that can pair with it (for some combination choice for the 2nd events) and give 
us solutions. We repeat this for every combination choice for every event. Neglecting ef- 
fects 2.- 7., Npair{c,i) = Nf.yt — 1 if c is the correct combination for event i. After including 
backgrounds and smearing, Npair{c,i) < N^vt — ^, but the correct combinations still have sta- 
tistically larger Npair{c,i) than the wrong combinations. The frequency with which various 
values of Npair{c, i) occur is shown as a function of Npair{c, i) in Fig. [61 

To enhance the likelihood that a particular choice of c, i corresponds to a correct solution, 
we cut on Npair{c,i). For the SPSla model point, if Npair{c,i) < 0.75 N^vt we discard the 
combination choice, c, for event i. If all possible c choices for event i fail this criterion, then 
we discard event i altogether (implying a smaller N^vt for the next analysis cycle). We then 
repeat the above procedure for the remaining events until no combinations can be removed. 
After this, for the example data sample, the number of events is reduced from 622 (424 
signal + 198 background) to 430 (322 signal + 108 background), and the average number 
of combinations per event changes from 11 to 4. 

Second, we increase the significance of the true solution by weighting each surviving 
pair of events by 1/n where n is the number of solutions for the given pair (using only 
the combination choices that have survived the previous cuts). This causes each pair (and 
therefore each event) to have equal weight in our histograms. Without this weighting, a pair 
with multiple solutions has more weight than a pair with a single solution, even though at 
most one solution would be correct for each pair. 

Finally, we exploit the fact that wrong solutions and backgrounds are much less likely to 
yield M^, Mx, My, and Mz values that are all simultaneously close to their true values. We 
plot the 1/n- weighted number of solutions as a function of the three mass differences (Fig. [7]). 
We define mass difference windows by 0.6x(peak height) and keep only those solutions for 
which all three mass differences fall within the mass difference windows. The surviving 
solutions are plotted (without the 1/n weighting) in Fig. [HI Compared with Fig. [5l the mass 
peaks are narrower, more symmetric and the fitted values are less biased. The fitted masses 
are {93.9, 140.3, 180.5, 559.2} GeV. Repeating the procedure for 20 data sets, we find 

rriN = 93.8 ± 3.9 GeV, mx = 138.4 ± 4.5 GeV, 

my = 178.7 ±4.6 GeV, = 559.5 ± 5.4 GeV , (24) 

to be compared to the input masses of Eq. f[22|) . Thus, the biases are reduced without 
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FIG. 7: SPSla, L = 300 fb~^ mass difference distributions. All effects incorporated, including 
backgrounds. 
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FIG. 8: Final mass distributions after the bias reduction procedure for the SPSla SUSY model 
and L = 300 fb~^. All effects incorporated, including backgrounds. 

significantly increasing the statistical errors. 

Thus, wc have shown that the masses can be measured with high precision for a few 
hundred events in the 4-fermion decay channel. In the case of the SPSla point, the number 
of events employed above corresponds to a high integrated luminosity, L ~ 300 fh~^. The 
reason that such a high luminosity is required in the case of the SPSla scenario is that the 
branching ratio for —>■ rr is 14 times that for — >■ /I// or — >■ ee. More generally, 
the integrated luminosity needed to get a few hundred events is highly dependent on the 
branching ratios for the various SUSY particle decays in the model. For example, if one takes 
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FIG. 9: Mass distributions for 50 events for SPSla. 

the SPSla masses but requires that X2 decays equally to the three lepton flavors instead, 
the same number of signal events as employed above can be obtained with just 10 fb~^ of 
data. 

Although the errors in the mass determinations depend upon the number of events, our 
method is quite robust in that we get decent mass determinations even with a small number 
of events. In Fig. [9], the mass distributions for 50 events are shown, with evident mass 
peaks. By repeating our procedure for multiple datasets of a given size, we obtain the errors 
as functions of the number of events. Fig. [10] shows the error for the Xi mass determination 
as a function of the number of signal+background events. Note that the central value for 
multiple data sets of the given size is quite insensitive to the data set size, but, of course, 
the possible deviation from this central value for any one data set increases as the data set 
size decreases. 




B. SUSY Point #1 

We have applied our method to other mass points to show its reliability. We quote here 
results for "point 7^1" defined in Ref. 6[] with the following masses: {85.3, 128.4, 246.6, 
431.1/438.6} GeV. For 100 fb~^ data, we have about 800 events (770 signal events) after 
the same pre-bias-reduction cuts. The resulting mass plot before performing bias reduction 
cuts is that given in Fig. [TTl From Fig. [11], we see that the mass peaks are very broad 
and we get more than 50 GeV biases if we use the positions of the maxima as the true 
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FIG. 10: Error bars for rn-jy as a functfon of the number of background+signal events, for SPSla. 
All effects and procedures included. 

mass values. We then repeat the same bias reduction procedure as for SPSla except that 
we employ a looser cut on Npair{c,i) than for the SPSla case, despite the fact that there 
are more signal events for Point #1. We require Npair{c,i) > 0.6 N^vt- The reason is that, 
unlike the SPSla case, the gluino mass in Point #1 (524 GeV) is significantly larger than 
the squark mass. Therefore the quark jet from gluino decay is often misidentified as the 
jet from squark decay, which reduces the chance to obtain solutions for a pair of events. ^ 
In practice, there is not a universal "best" cut on Npair{c,i): a more stringent cut leads to 
smaller biases but larger statistical uncertainties. After the bias reduction procedure using 
Npair{c,i) > 0.6 Nevt we are left with 560 events (550 signal events). The mass distributions 
are shown in Fig. [121 They are much narrower and the biases are considerably reduced. 
After following the bias reduction procedure and using 20 data samples to estimate the 
errors, we obtain = 82.8 ± 3.2 GeV, mx = 127.9 ± 3.0 GeV, my = 245.7 ± 3.4 GeV, 
niz = 436.4 ± 5.4 GeV. The central values are in quite close agreement with the input 
masses except for rriN- which comes out a bit low. 



^ It is possible to improve the results by considering all high pr jets as candidates for the quarks from 
squark decays, instead of simply choosing the two highest px jets. 
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FIG. 11: Final mass distributions before the bias reduction procedure for the point #1 SUSY 
model and L = 100 fb^^. All effects incorporated, including backgrounds. 
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FIG. 12: Final mass distributions after the bias reduction procedure for the point #1 SUSY model 
and L = 100 fb~^. All effects incorporated, including backgrounds. 

C. Comments and Comparisons 



We emphasize that the remaining biases in the above mass determinations can be removed 
by finding those input masses that yield the observed output masses after processing Monte 
Carlo generated data through our procedures. In this way, very accurate central mass values 
are obtained with the indicated statistical errors. 

The above results for the A^, Y and X masses for the SPSla point and point #1 can be 
compared to those obtained following the very different procedure of Ref. 6|. There, only 
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the X ^ Y ^ N parts of the two decay chains were employed and we used only 4/i events. 
For the SPSla model point we obtained uin = 98 ± 9 GeV, my = 187 ± 10 GeV, and nix = 
151 ± 10 GeV. And, for point #1 we found m^v = 86.2 ± 4.3 GeV, nix = 130.4 ± 4.3 GeV 
and my = 252.2 ± 4.3 GeV. Including the 4e and 2/i2e channels will reduce the indicated 
errors by a factor of ~ 2. The procedure of [Q] can thus be used to verify the results for m^v, 
mx and my from the present procedure and possibly the two can be combined to obtain 
smaller errors than from either one, with mz determined by the procedure of this paper. 

We also compare the results for SPSla with those given in Ref. {t] where exactly the 
same procedure and cuts are applied to the same model point. The difference is that we 
used ATLFAST for the detector simulation in Ref. [3] while we have switched to PGS in the 
current paper. The PGS simulation has more stringent lepton isolation cuts and therefore 
we obtain fewer events in the present analysis (620 vs 1050). In Ref. p], we obtained 

m^ = 76.7 ± 1.4 GeV, m^ = 135.4 ± 1.5 GeV, 
my = 182.2 ± 1.8 GeV, mz = 564.4 ± 2.5 GeV. 

before the bias reduction procedure and 

m7v = 94.1 ±2.8 GeV, mx = 138.8 ± 2.8 GeV, 
my = 179.0 ± 3.0 GeV, mz = 561.5 ± 4.1 GeV. 



after. Comparing the above numbers with those in Eqs. fl23|) and fl24|) . we see that the 
masses obtained using PGS simulation have larger statistical errors, in accord with the 
smaller number of events. On the other hand, the central values agree well, indicating that 
the bias reduction procedure affects the mass peaks in a nearly model-independent manner. 
We view this as evidence of the robustness of our method. 



D. Removing Biases Using a Dilepton Edge Cut 

As we have discussed, the primary source of biases is the detector smearing of wrong 
solutions, especially those associated with wrong combinations. It will be possible to effi- 
ciently eliminate many of these wrong solutions if there is a significant structure associated 
with correct solutions in one or more distributions constructed from the visible particles' 
momenta. In the SUSY examples we consider here, such a structure is especially apparent 
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FIG. 13: We plot the number of events as a function of for 600 SPSla events in 5 GeV 

bins (after PGS smearing and general cuts, but before the bias reduction procedure). Only events 
containing two muons and two electrons with opposite charges are used to avoid ambiguity, each 
of which conntributes two entries to the histogram. The edge at 80 GeV is apparent. 

in the distribution of m£+£-, where i = e,fi (same flavor pairs only). The advantage of 
using only leptons is the much better resolution for the lepton momentum measurements. 
Ignoring resolution smearing, kinematics predicts that correct combinations should have 



/2 2 \ / 2 2 \ 

/ edge \2 _ l^yp ~ ^Xo)\^XO ~ ^No) 



m 



(25) 



xo 



where myo-, ^xo and niNo are the input masses. Note that there are many more dilepton 
events than four lepton events since dileptons require only a single decay chain of Fig. [H The 
plot of m£+£- values for all solutions coming from 600 SPSla events (after PGS smearing 
and general cuts, but before applying the bias reduction procedure) is shown in Fig. [131 The 
edge at the predicted value of 80 GeV is apparent and its location can be determined quite 
accurately from the data. 

Before employing the bias reduction procedure, we apply a cut on the my, mx, tun values 
obtained for a given solution of 



m 



< 20 GeV, 



(26) 



where we have purposely employed a rather loose cut so as to not lose statistics and to 
take into account smearing of the input X, Y, N masses that will be present even for a 
correct combination, as well as the small error associated with determining the edge location 
experimentally. We then use the same bias reduction procedure as discussed earlier using a 
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sequence of choices for the cut fcut defined by retaining only combinations with Npair{c, i) > 
fcutNevt- In Table [H 





with dilepton edge cut 


without dilepton edge cut 


fcut 


0.60 


0.65 


0.70 


0.75 


0.60 


0.65 


0.70 


0.75 


TUN (GeV) 


93.0 ±3.7 


96.1 ±3.9 


97.5 ±4.3 


97.9 ±4.9 


85.6 ±2.3 


88.1 ±3.5 


90.7 ±3.8 


93.8 ±3.9 


mx (GeV) 


138.9 ±3.9 


141.4 ±4.6 


143.7 ±4.6 


144.3 ±4.0 


131.5 ±2.7 


133.9 ±3.6 


135.9 ±4.3 


138.4 ±4.5 


my (GeV) 


176.5 ±3.8 


178.8 ±4.6 


180.8 ±5.1 


181.5 ±5.3 


172.8 ±2.8 


174.8 ±3.8 


176.6 ±4.4 


178.7 ±4.6 


mz (GeV) 


557.8 ±4.4 


559.9 ±4.5 


563.2 ±5.0 


565.6 ±6.2 


555.8 ±5.2 


557.2 ±5.5 


557.8 ±5.1 


559.5 ±5.4 



TABLE I: Peak locations for various values of fcut with and without the dilepton edge cut. Errors 
were determined using 20 distinct data sets. 

We clearly observe that the dilepton edge cut has greatly reduced the bias in comparison 
to results obtained without the dilepton edge cut. Further, for the larger values of fcut, the 
mass peak locations are not biased at all (within statistics) in comparison to the input masses 
of Eq. ( 122|) . This occurs because many of the wrong solutions have been elliminated. For 
example, after the dilepton edge cut and after employing the bias-reduction procedure using 
fcut = 0.75, about 160 events are retained on average and the average number of solutions for 
the remaining pairs formed from these surviving events is only about 1.2. In other words, the 
dilepton edge cut is highly effective in removing wrong combinations. Errors for the peak 
mass values are, of course, slightly larger when a dilepton edge cut is imposed, implying 
that ultimately the best mass determinations may be those obtained using Monte Carlo 
determination of the bias corrections that should be applied to mass peak values obtained 
without the dilepton edge cut. Nonetheless, doing the analysis with a dilepton edge cut will 
provide a very important cross check of the bias determination. 

As a final note, we observe that in the case of SPSla there are also incorrect solutions 
coming from chains containing a pair of leptonically decaying r's. Many, but not all of these 
wrong solutions are also be eliminated by the dilepton edge cut. The remaining background 
events contain mostly those events for which one chain has £ = e, /i while the other has a 
pair of leptonically decaying r's, since sometimes such events will give a solution with nearly 
correct mass values. 
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FIG. 14: Dilepton invariant mass distribution for 600 SPSla events together with 600 Point 1 
events (after PGS smearing and general cuts, but before the bias reduction procedure). Only 
events containing two muons and two electrons of opposite charge are used to avoid ambiguity. 
The two edges at 80 GeV and 157 GeV correspond to SPSla and Point #1 respectively. 

E. More on Backgrounds 

Because the SM background can be efficiently reduced by applying a large missing px 
cut, the most difficult backgrounds usually come from other SUSY processes that contain 
the same final state particles. In the above examples, we have already encountered such 
backgrounds. In the SPSla case, the backgrounds are dominated by events that contain 
leptonically decaying r's. For SUSY point #1, although most events are originally signal 
events, in many cases the jets from squark decays arc not correctly identified, in which case 
these events should be viewed as background events. In both examples, the background 
events are closely related to the signal events and therefore also carry some information 
about the masses. It is also interesting to study the effects of background events of a 
completely different origin, and test the stability of our mass determination method. 

In order to explore the issues that arise, we will perform an analysis in which we consider 
the SPSla events as signal events, fixing the number of events to 600 (including the intrinsic 
background of SPSla). For a possible SUSY background to the SPSla events we employ 
events of the above SUSY Point #1 as "background" events, varying the ratio of Point 
#1 events with respect to the SPSla events. Since this is only for illustration, we are not 
concerned with how this could happen in a specific SUSY model. The existence of two 
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different type of events is immediately seen from the dilepton invariant mass distribution 
(Fig. [T^ . where two different edges are evident. The position of the edges are given by 
Eq. (12^ . In this case, the two signals give two different myo, mxo and ttinq input mass 
sets. Again, there are many more dilepton events than four lepton events since a dilepton 
only requires a single decay chain of Fig. [T] From Fig. [T3] we see that the position of the 
edge associated with the SPSla signal can be determined quite precisely even when the 
background to signal ratio is of order one. Consequently, one can try to combine the edge 
location measurement with information from double chain events (see also Ref. 2l|). 



First, we repeat our fitting procedure on the mixed events without using the dilepton edge 
information. Since there are more background events, we cannot use the fixed Npair{c,i) > 
0.75 Ne^t cut as before. This is because N^vt now refers to the total number of events from 
both SPSla and SUSY Point ^1 so that such a cut would amount to a much stronger 
fcut value for the SPSla signal on its own. Instead, we choose the Npair{c,i) cut so that 
60% of all the events are left after the bias reduction procedure. The corresponding value 
of fcut varies according to the amount of SUSY Point #1 background included. For 600 
SPSla events combined with 600 Point ^^1 background events the 60% survival fraction 
corresponds to using fcut ~ 0.58 in the bias-reduction procedure. For the SPSla signal 
alone, the corresponding fcut value is somewhat larger. The measured m^r is shown in 
Fig. [15] as a function of the number of background events (after PGS smearing, but before 
the bias reduction procedure). From Fig. [15], we see that the mass determination is not 
accurate when the background/signal ratio is high. However, as long as the number of 
background events after general cuts is less than about half the number of signal events, the 
bias reduction procedure is effective in removing background events while retaining signal 
events, and the mass determination is quite good. Of course, in practice we will not know 
a priori what the number of background events is relative to the number of signal events 
and therefore we would need additional input in order to know if the fitted masses from 
mass peaks are reliable estimates for the true masses. In the present case, the dilepton mass 
plot of Fig. [14] would clearly have indicated the presence of two different classes of events 
and we would therefore know that it would help to use an additional cut to reduce the 
"background" class. We again employ the simple cut of Eq. (1261) . where in order to isolate 
the SPSla component of the combined events we would employ "^^+^- = 80 GeV. Again 
note that this is a loose cut that does not require a precise knowledge of the edge position 
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FIG. 15: The measured m^v as a function of the number of background events. The number of 
signal (SPSla) events is fixed to 600. In both cases, these are the event numbers after general cuts 
but before the bias reduction procedure. 

from the single chain events. Applying the cut in Eq. (126|) on datasets with 600 signal + 
600 background events (after general cuts, but before the bias reduction procedure) and 
repeating the fit procedure, we obtain the masses 



where errors were determined using 20 distinct data sets. The same fcut = 0.58 was used 
as in the case without dilepton edge cut. With the dilepton edge cut, we are left with 
averagely 363 SPSla events and 8 SUSY Point 7^1 events in the mass distributions used to 
get the mass peak locations of Eq. fl27|) . Thus, we effectively isolated the signal of interest 
by employing the dilepton edge cut. We have also obtained central mass values that have 
almost no bias relative to the input SPSla masses. This is because the fcut for the SPSla 
signal alone that yields the number of events (~ 360) after bias reduction is close to 0.65, 
which according to Table H] should give a nearly unbiased mass determination. 

V. COMPARISON OF SUSY AND UED 

In this section, we address the question of whether or not the mass determinations (and 
the accuracy thereof) are sensitive to the model employed by comparing results for the SPSla 



TUN = 96.2 ± 3.6 GeV, nix = 141.3 ± 4.1 GeV, 
my = 178.4 ± 4.1 GeV, mz = 558.5 ± 4.6 GeV 



(27) 
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FIG. 16: Final mass distributions for signal events only after employing the bias reduction proce- 
dure for the SPSla mass choices in the context of the SUSY model and using 400 signal events 
after PGS smearing and general cuts but before bias reduction. All combinations and solutions 
are included, but backgrounds associated with the SUSY model are not included. 

point to a UED model chosen so that the masses are exactly the same as for SPSla (the 
corresponding decay chain in UED is: KK-quark — )■ KK-Z — )■ KK-lepton — )■ KK-photon). We 
have also adjusted the squarks/KK-quarks of different flavors to have the same mass (564.8 
GeV, as for ul in SPSla) and chosen squark/KK-quark pair production as the only process 
{i.e., no gluinos/KK-gluons). The finite widths of the involved particles are also turned 
off. Both SUSY and UED events are simulated with Herwig++ [39], with spin correlations 
included and confirmed by comparing with Ref. 40[. In Figs. [16] and [T71 for SUSY and 
UED respectively, we plot the mass distributions after employing identical smearing (PGS), 
general cuts and the bias reduction procedure (using Npair{c, i) > 0.75 Nevt)- The size of the 
event samples for SUSY and UED are set by requiring that both samples contain 400 events 
after PGS smearing and general cuts, but before the bias reduction procedure. Visually, it 
is clear that the peaks are in very similar locations. 

After employing our standard fitting techniques, we obtain masses of 



= 90.3 ± 3.0 GeV, = 135.4 ± 3.1 GeV, 

my = 176.0 ± 2.9 GeV, mz = 551.2 ± 5.2 GeV. 



for the actual SUSY SPSla model and masses of 



niN = 89.5 ± 3.5 GeV, mx = 134.8 ± 3.5 GeV, 
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FIG. 17: Final mass distributions for signal events only after the bias reduction procedure for the 
SPSla mass choices, but in the case of the UED model. The UED event sample is scaled so that 
there are 400 signal events after PGS smearing and general cuts but before bias reduction. All 
combinatorics and solutions are included, but backgrounds associated with the UED model are not 
included. 

my = 176.2 ± 3.7 GeV, mz = 548.7 ± 5.8 GeV. 

for the UED model with SPSla masses. We note that the measured masses still have some 
biases of order a few GeV (for the three smaller masses) up to more than 10 GeV. But the 
biases are very similar for the two models with different spins, indicating that the major part 
of these biases can be removed by comparing with Monte Carlo even before the underlying 
model is determined. Using the technique of repeating the procedure 20 times, we find that 
the SUSY model statistical errors in this case are somewhat smaller than quoted earlier in 
Eqs. ( 124|) because we have not included any background from other SUSY processes. The 
biases are also slightly different from Eqs. fl24|) where the backgrounds are included. The 
UED model errors would also be increased by including backgrounds coming from other UED 
processes. In other words, the masses are very well determined (to within a few GeV) by our 
purely kinematic procedures, but the errors are mildly model-dependent because of variation 
in the nature and magnitude of the new-physics-model backgrounds. The contribution from 
the backgrounds can often be inferred from real data and subtracted. For example, the 
SPSla background is dominated by the chain decays that yield stau's, whose existence and 
production rate can be determined from hadronically decayed tau's. 
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VI. DISCUSSION AND CONCLUSIONS 



One important question is whether it is better to use one-chain or two-chain techniques. 
Our point of view is that one should use all available kinematic information, regardless of 
whether it is from one-chain or two-chain events. On the one hand, due to the fact that 
not all events contain two identical chains, one often obtains more one-chain events of a 
certain type than events with two identical decay chains of that type. However, if one 
considers only one chain at a time, information, in particular that related to the measured 
missing transverse momentum, is always lost. The consequence is that either one cannot 
solve directly for all involved masses for a given length of decay chain, or one must employ 
longer decay chains, in which case the method becomes very complicated. For example, in 
the one-chain case one needs to employ decay chains with four visible particles (vs. three 
visible particles in the two-chain case) and, in addition, one needs to combine five events to 
obtain discrete solutions for the unknown masses. Even assuming that such events do exist, 
there are more wrong combinations and wrong solutions than the two-chain case studied in 
this paper. Further, the existence of a certain type of decay chain implies that there are 
always events with two identical such decay chains. Events with two identical decay chains 
always provide more information for the masses of the particles in the decay chains. The 
challenge of two-chain techniques is that one needs to identify those events in which there 
are indeed two identical chain decays. Ideally, one would divide the observed events into 
different channels according to their event topologies (chain type 1 -|- chain type 1, chain 
type 1 + chain type 2, chain type 2 + chain type 2, . . . ), apply methods appropriate to each 
topology and, in order not to lose statistics, combine all these channels in the analysis. For 
this reason, it is very important to extend the studies in Ref. {g] and the current paper to 
other event topologies. 

The importance of using one-chain decay information is illustrated in the SPSla case. 
Since there are many more dilepton events than 4-lepton events the dilepton edge given 
in Eq. ( 125|) can be measured very precisely. If available, one should certainly incorporate 
this measurement into the two-chain techniques to better determine the masses. This kind 
of "hybrid" approach has been studied here and in {20, 
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231 . As summarized below. 



adding the dilepton edge information improves the two-chain mass determinations obtained 
following the basic procedure developed in this paper. 
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In our procedure, we have applied a set of bias reduction methods. In particular, to reduce 
the number of combinations, we have utilized the fact that correct combinations can pair 
with relatively more events than wrong combinations. Alternatively, one may try to reduce 
the number of wrong combinations before doing the event reconstruction. For example, one 
could try to group objects into two hemispheres [41] and assume that the objects in the 
same hemisphere come from the same decay chain. However, this only works well when the 
initial particles are substantially boosted, while the squarks in our case are produced mostly 
close to the threshold without large boosts. For small boosts, the quark and X2 from squark 
decay actually belong to two opposite hemispheres instead of the same one. The directions 
of the subsequent X2 ^ decay products are even more random. We have applied the 
hemisphere method on a set of ideal events from squark pair production, each containing 2 
quarks and 4 leptons according to the decay chain in Fig. [H Even without any complications 
from extra jets, experimental smearing and so forth, only about 12% of the events have the 
decay chains correctly identified (this does not account for the ambiguity of the two leptons 
in the same decay chain). 

One could explore the effectiveness of imposing a cut which accepts only events with 
substantial thrust or small circularity before separating each event into two hemispheres. It 
is not clear to us that the gain from decreasing the combinatorics problem would outweigh 
the reduced statistics associated with the fact that such a cut would remove a large fraction 
of the available events. 

However, we have shown that when a two-particle mass edge (such as the dilepton mass 
edge in the examples we have considered) can be identified, it is very useful to impose a 
cut whereby only solutions that give a dilepton mass within roughly ±20 GeV of the edge 
location are retained. By applying this cut before proceeding with the rest of our analysis 
procedure the mass peak biases are essentially eliminated and the errors on the central mass 
values are very similar. In addition, we have shown that this cut is capable of essentially 
eliminating contamination of the mass determinations for the SUSY signal of interest by 
events coming from some other SUSY signal that does not share the same dilepton edge 
location. Presumably any other recognizable kinematic edge could be exploited in similar 
fashion, but dilepton mass edges will typically be least impacted by detector momentum 
smearing. 

In conclusion, we have proposed a kinematic technique for mass determination in events 
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with two invisible dark matter particles. The technique seeks constraints on the mass space 
from measured momenta. In Sec. [Tll we have given general constraint counting and discuss 
the corresponding strategies for both the single decay chain case and the double decay chain 
case. In the former, one only uses the information from one of the two decay chains in 
each event; in the latter, one uses information associated with both decay chains. The 
constraints include the mass-shell constraints for the dark matter particle as well as all 
intermediate particles. In the double decay chain case, we obtain extra constraints from the 
measured missing transverse momenta. In both cases, more constraints are available when 
the decay chain is longer. In certain instances (this includes the single decay chain case with 
3 visible particles per chain and the double identical decay chain case with 2 visible particles 
per chain), we obtain discrete solutions for the missing momenta by using trial masses for 
the unknown particles. Requiring that we obtain physical solutions for the momenta, the 
consistent mass region becomes more restricted when more events are included. The actual 
masses are obtained by studying the kinematic distributions of the consistent region. This 
is the strategy we adopted in Refs. jol, |23| . When the decay chains present in each event are 
longer, it is possible to obtain discrete solutions for the momenta (and therefore the masses) 
by combining the constraints from a few different events, without assuming any trial masses. 
This occurs when there are 4 or more visible particles per decay chain for the single chain 
case, or 3 or more visible particles per decay chain for the double chain case. In this article, 
we have focused on the double identical chain case with the number of visible particles per 
chain fixed to 3. 

In our case, the constraints can be solved for discrete solutions of the unknown masses 
when two events are combined. However, because the system of equations contains quadratic 
equations, wrong solutions are introduced. Nonetheless, if the visible momenta could be 
measured without errors, it would take only three events to obtain the correct masses. This is 
because the wrong solutions would be different for the three different possible event pairings 
whereas the correct solution remains the same for every event pair. This remains true even 
when wrong combinations are included. However, in practice the non-zero experimental 
resolutions imply that the correct solution distribution becomes smeared which then overlaps 
with the distribution coming from wrong solutions, wrong combinations and background 
events. Despite this, we have shown that when two of the visible particles in a decay chain 
are leptons, we obtain good precision (~ a few GeV) if a few hundred events are available. 
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We have developed methods to reduce the number of wrong combinations and backgrounds. 
The resulting mass solution distributions are clearly peaked around the input masses with 
small systematic errors which can be eliminated either by comparing with Monte Carlo 
distributions around the estimated masses or by imposing an initial dilepton edge mass cut 
on the accepted solutions. An important assumption for doing this comparison is that the 
distributions are only sensitive to the masses instead of the underlying theories. We have 
shown that this is indeed the case by comparing two distinct theories with different spin 
structure, MSSM and UED. We have set the spectra of the two models to be the same 
and examined the mass distributions, which show little difference. Correspondingly, the 
systematic errors introduced by model dependence are much smaller than the statistical 
errors. 

Finally, we comment on possible improvements of our method. Given theprecision of the 
purely kinematic results and availability of spin determination techniques 42(| it would cer- 
tainly be possible to figure out the underlying theory/spins and then apply model-dependent 
techniques to refine the mass determinations. For example, one could adopt a likelihood 
method similar to the one used in the top mass measurement with dilepton events at Teva- 
tron [43]. In this method, a probability density, as a function of all unknown masses, is 
defined for each event by convoluting matrix elements and detector resolution functions. 
One then obtains the joint probability by taking the product of the probability densities 
from all events. The best-fit masses are given by the values that maximize the joint proba- 
bility. Compared with top quark mass measurement, the event topology considered in this 
paper is more complicated since that more final state particles are involved. One also needs 
to scan a four dimensional instead of one dimensional mass space to minimize the probability 
function. Therefore, another benefit of the purely kinematic method is that it significantly 
simplify the computation by reducing the candidate mass space to a very small region. 

We could also consider simplified likelihood methods using the detector resolution func- 
tions, but without the knowledge of the matrix element. One such method is based on the 
same constraints, Eqs. (E]) through ( ITOi) . which can be viewed as a generalization of the 
method discussed in this paper: for each event, we can eliminate the 4-momenta of the two 
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missing particles by using Eqs. (l6|)- (fT0|) . We are then left with two equations in the form 

f{mN,mx,mY,mz]Pvis) = 0, (28) 
g{mN,mx,mY,mz;Pvis) = 0, (29) 

where / and g are functions of the masses and all visible momenta, Pms- If we require that 
the equalities hold exactly, then each event defines a 2-dimensional surface in a 4-dimensional 
mass space. The 2-dimensional surfaces of two different events intersect at discrete points, 
which are nothing but the mass solutions. This is equivalent to saying that in the mass 
space, for each event we assign a non-zero probability density for points on the surface and a 
zero probability density for points off the surface. Then by combining two events, the joint 
probability is non-zero only at discrete points. Obviously, a more sophisticated method is 
to assign a maximum probability density for points on the surface, and smaller but non-zero 
probability densities for points away from the surface. This can be done by calculating the 
distribution for each point in the mass space, 

= + f,, (30) 



- E ( ) . (31) 



where 

and a similar formula holds for cr^. Here, we have assumed Gaussian distributions for 
the invisible momenta Pms with errors given by cip^.^, and Eq. f l3T|) can be complicated by 
correlations among the visible momenta. We then sum the over all events and obtain the 
masses when the total is minimized. In practice, / and g are very complicated functions 
of the visible momenta, which may result in large roundoff errors. As in the matrix element 
method, one also needs to efficiently minimize the over a 4-dimensional space. Therefore, 
it deserves further studies to determine whether this approach can improve the precision of 
the mass measurement over the simple method discussed in this paper. 

In closing, we note that the program for determining the Tnz,Y,z,N solutions 
as a function of input visible momenta for two-chain events is available at 



http: //particle. physics. ucdavis.edu/hefti/projects/ 1 (choose WIMPMASS). Stand-alone pro- 
grams that implement the methods of Refs. 6|, 123] are also available at this same website. 



34 



Acknowledgments 



This work was supported in part by the U.S. Department of Energy grant No. DE-FG02- 
91ER40674. 



Appendix A: Solving the constraint equations 

In this Appendix we describe in detail our procedure to solve the system of the kinematic 
constraint equations 

pi = pl = ql = ql (Ai) 

2pi -^3+^3 = 2p2 ■P4+pI = 2gi ■q3 + ql = 2g2 ■ q^ + ql, (A2) 
2(pi + Ps) ■P5+pI = 2(p2 + P4) ■P6+pI = 

= 2(gi + ^3) • + ql = 2(g2 + ^4) ■ ge + ql, (A3) 
2(pi +P3+ Ps) ■P7+P7 = 2(p2 +P4 + Pe) ■P8+pI = 

= 2{qi + gs + gs) ■ g? + g? = 2(g2 + g4 + ge) ■ gs + gs, (A4) 

Pl+P2= Pmiss^ Pl+P2= Pmiss^ (A5) 

g? + g2=gmis., gf + g| = gLss- (A6) 

As mentioned earlier, it is straightforward to numerically solve the above equations using 
commercial software such as Mathematica, but the speed is intolerably low. Instead, we 
solve the equation system using a programming language such as C++. The idea is to reduce 
the system to a univariate polynomial equation whose coefficients are (fixed) functions of 
the original visible momenta. Note that it is convenient to obtain the coefficient functions 
with the assitance of Mathematica. After that, the functions are hard coded in the C++ 
program and Mathematica is no longer needed. The univariate equation can then be solved 
numerically using any available polynomial solver. The Mathematica notebook and the C++ 
code are available from Ref. {44] or any of the authors. The key Mathematica operations are 
also described in Appendix [B] for reference. The method can be easily generalized to solve 
other polynomial equations efficiently. 

It is straightforward to eliminate 13 variables using the 13 linear equations and obtain 3 
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quadratic equations with 3 variables. Generically, 3 quadratic equations can be written as 

+ a^zy + a-jzx + + a^ip' + ai^yx + a^y + 0.2^^ + oi^: + ao = 0, (A7) 
2;?/ + 672;x + h^z + 65?/^ + h^yx + 63?/ + 622;^ + h^x + 60 = 0, (A8) 
2:0; + cqz + C5?/^ + Ci^yx + C3?/ + C2X^ + c\x + Co = 0. (A9) 

where x,y,z are variables, and ai,bi,Ci are coefficients as functions of the original visible 
momenta. We have ordered the left-hand terms lexicographically in the order z > y > x. 
We will eliminate variables also in this order and eventually obtain a univariate equation in 
X. In our implementation, we choose x, y, z to be Pi,P2, Qi or some permutation thereof. In 
fact, we solve several times for several different permutations just to make sure we do not 
miss any solutions. 

First, by calculating (]A7p x y — (lASP x z, we cancel the term z'^y and obtain a polynomial 
equation with leading term —h-jxz"^. We repeatedly use flA7p . flASp . flA9P to reduce the 
polynomial, i.e., to eliminate the leading term of the polynomial. For example, we eliminate 
the —bjxz'^ term by subtracting —bj x ( 1A9|) . The next leading term is oc z"^, which again 
can be eliminated by subtracting from it Eq. ( 1A7I) with appropriate coefficient. Repeating 
this procedure until it cannot be reduced further, we obtain an equation in the form 

z + y^ + y^x + y"^ + yx^ + yx + y + x'^ + x"^ + X + 1 = 0, (AlO) 

where we have omitted all coefficients. Similarly, reducing ( 1A7P x x — ( ]A8[) x 2; , we obtain 
another equation also in the form of flAlOp . Canceling the leading term, we eliminate the 
variable z and obtain a cubic equation in {y,x). Reducing (lASP x x — (1A9P x y, we obtain 
another cubic equation in {y,x). 

We are then left with 2 cubic equations in 2 variables, 

y"^ + dgy'^x + djy"^ + d^yx^ + d^yx + d^y + d^x^ + d2X^ + dix + c/q = 0, 

y X + e-jy + e^yx + e^yx + e^y + e^x + e2X + Cix + Cq = 0, (All) 

where the coefficients di and Cj are derived from Oj, bi and q in Eq. (IA7D ( 1A8I) ( 1A9I) . This 



system can be solved using the following method. The resultant [29| of two univariate 
polynomials is defined as follows. Given a polynomial (note that the and bi below are not 
those given in Eq. flA7p and flASP ) 

P{x) = anx"" + a„_ix""^ + ... + aix + ao , (A12) 
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of degree n with roots 



1 , . . . , n and a polynomial 



Q{x) = hmx"" + hm-ix"^-' + ... + hix + 60 



(A13) 



of degree m with roots j = l,...,m, the resultant p{P,Q), also denoted R{P,Q) and 
also called the eliminant, is defined by 



piP,Q) = a^blllllia,-P,) 



(A14) 



30|. 



The resultant is also given by the determinant of the corresponding Sylvester matrix 
The Sylvester matrix associated to polynomials P and Q is the (n + m) x {n + m) matrix 
obtained as follows: 



1. the first row is: 



(^ttn ■ ■ ■ ai ao ■ ■ ■ 



(A15) 



2. the second row is the first row, shifted one column to the right; the first element of 
the row is zero. 

3. the following (m — 2) rows are obtained the same way, still filling the first column with 
a zero. 



4. the (m + l)-th row is: 



bm. b, 



m "m— 1 



61 60 ■ ■ ■ . 



(A16) 



5. the following rows are obtained the same way as before. 

Taking the left-hand side polynomials of Eq. (lAlip as P and Q respectively, we have n = 3 
and m = 2 and the Sylvester matrix in this case is 

^ 03 02 Oi ^ 

03 02 fli flo 

S= 62 &i &o , (A17) 

62 bi bo 

\ b2 b, bo J 



where 



03 = 1, a2 = rfgx + dj, ai = d^x"^ + d^x + ^4, ao = d-^x^ + d2x'^ + dix + do, 

b2 = x + 67, bi = cqx^ + 650; + 64, 60 = 630;^ + 620;^ + eix + cq. (A18) 
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Therefore, p(P, Q) = detS is a 9th order polynomial in x. By Eq. (lA14p . if p(P, Q) = 0, the 



two equations in flAll|) have at least one common root. In other words, for each root Xj of 



detS* = 0, we can find a i/i such that {xi,yi) is a solution of Eqs. (lAlip . Thus, the problem 
has been reduced to the solution of a 9th order polynomial equation with real coefficients. 
We solve it numerically using algorithm TOMS/493 Note that one of the roots is fake. 
It can easily be identified and eliminated by substituting back each of the solutions into the 
original equations to identify the root that is not actually a solution of the original equations. 

Appendix B: Mathematica file 

Many of the operations described in Appendix |A] can be conveniently done in Mathemat- 
ica, which we describe below. 

First, suppose we have obtained the equations flA7p . flASp and flA9p . As mentioned, the 
coefficients Oj, bi and q are functions of the visible momenta. These functions are usually 
complicated, therefore it is always desirable to use intermediate parameters such as Oj, bi 
and Cj, and normalize the coefficients of the leading terms to 1. To obtain Eqs. (lAlip . we 
use the Mathematica function PolynomialReduce: 

PI = aO + al X + a2 x'^ + a3 y + aAy X + a5 y'^ + a6 z + a7 z X + a8 z y + z"^ ; 
P2 = b0 + blx + b2x'^ + b3y + bAyx + b5y^ + b6z + b7zx + zy; 
P3 = cO + cl X + c2 x^ + c3y + cAy X + cby"^ + c6 z + z x; 
PA = Expand[PolynomialReduce[Ply - P2z, {PI, P2, P3}, {z, y, x}][[2]]] 
P5 = Expand[PolynomialReduce[Plx - P3z, {PI, P2, P3}, {z, y, x}][[2]]] 
P6 = Expand[PolynomialReduce[P2x - P3y, {PI, P2, P3}, {z, y, x}][[2]]] 

where we obtain PA and P5 as two polynomials in the form of Eq. flAlOp and P6 a cubic 
polynomial in x and y. Then it is straightforward to cancel the variable z and cast the 
remaining two polynomials in the form of Eq. (lAlip . The resultant is obtained from 



P7 = d0 + dlx + d2x^ + d3x^ + dAy + d5yx + d6yx^ + d7y^ + d8y^x + y^; 
P8 = eO + elx + e2x +e3x + eAy + eby x + e6y x +e7y +y x; 



resultant = Resultant[P7 , P8, y]; 
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where resultant is the 9th order polynomial in x we are seeking. After obtaining the roots 
of the resultant numerically, y and z can be uniquely determined as follows: defining the 
polynomial P9 by 

P9 = Expand[PolynomialReduce[P7 X - P8y, {P7, P8}, {y, x}][[2]]]; 

we see that P9 can be written in the form 

P9 = /0 + /lx + /2x^ + /3x^ + /4x'^ + /5y + /6yx + /7yx^ + /8yx^ + y2; 

It turns out that the polynomial PIO defined below is linear in y: 

PIO = Expand[PolynomialReduce[P9x - P8, {P8, P9}, {y, 2;}][[2]]]; 

We then obtain y from the equation PIO = by substituting in the solutions for x. Once 
this is done, we can obtain z from the equation PA = since it is linear in z. One of the 9 
solutions is a fake solution to the original equation system, which can be easily identified by 
finding the solution which does not actually solve the system of equations when substituted 
back into the system of equations. 
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